www.gusucode.com > A benchmark software for MSPC工具箱matlab程序 > A benchmark software for MSPC/qFactor.m

    
function [sys,x0,str,ts] = qFactor(t,x,u,flag)
%qFactor_S_Function: Calculates q factor from composition z and temperature T

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes();

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);

  %%%%%%%%%%%%%%%%%%%
  % Unhandled flags %
  %%%%%%%%%%%%%%%%%%%
  case { 2, 4, 9 },
    sys = [];

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end
% end csfunc

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts]=mdlInitializeSizes()

sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = [];


%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim
%    state
% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)

sys = [];

% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)

%Data
%********************************
VLEdata = load('VLEdata.mat');
xA = VLEdata.xA;
Tbubble = VLEdata.Tbubble;
Tdew = VLEdata.Tdew;

M_methanol = VLEdata.M_methanol;
M_ethanol = VLEdata.M_ethanol;
CpL_methanol = VLEdata.CpL_methanol;
CpL_ethanol = VLEdata.CpL_ethanol;
CpV_methanol = VLEdata.CpV_methanol;
CpV_ethanol = VLEdata.CpV_ethanol;
Tn_methanol = VLEdata.Tn_methanol;
Tn_ethanol = VLEdata.Tn_ethanol;
Tc_methanol = VLEdata.Tc_methanol;
Tc_ethanol = VLEdata.Tc_ethanol;
Hv_methanol = VLEdata.Hv_methanol;
Hv_ethanol = VLEdata.Hv_ethanol;
%********************************


%Calculate Tbubble and Tdew from zF value:
z = u(1);
T = u(2);
Tbp = interp1(xA,Tbubble,z);
Tdp = interp1(xA,Tdew,z);

%Possible cases:
% Feed condition          q       Case
%---------------        -----   ------------
% Subcooled liquid       >1     (1) TF < Tb
% Saturated liquid        1     (2) TF = Tb
% Vapor-liquid mixture	0<q<1	(3) Tb < TF < Td
% Saturated vapor         0     (4) TF = Td
% Superheated vapor      <0     (5) TF > Td

if T<Tbp;
    %Subcooled liquid
    CpL = CpL_methanol*z+CpL_ethanol*(1-z);
    %MF = M_methanol*z+M_ethanol*(1-z);
    HL = CpL*(Tbp-T);
    %Corrección por temperatura:
    Hv_methanol = Hv_methanol*1000*((Tc_methanol-(Tbp+273.15))/(Tc_methanol-Tn_methanol))^0.38;
    Hv_ethanol = Hv_ethanol*1000*((Tc_ethanol-(Tbp+273.15))/(Tc_ethanol-Tn_ethanol))^0.38;
    HG = z*(CpL_methanol*(Tdp-T)+Hv_methanol)+(1-z)*(CpL_ethanol*(Tdp-T)+Hv_ethanol);
    q = HG/(HG-HL);
elseif T==Tbp;
    %Saturated liquid
    q = 1;
elseif Tbp<T && T<Tdp
    %Vapor-liquid mixture
    
elseif T==Tdp;
    %Saturated vapor
    q = 0;
elseif T>Tdp
    %Superheated vapor
    CpV = CpV_methanol*z+CpV_ethanol*(1-z);
    %Corrección por temperatura:
    Hv_methanol = Hv_methanol*((Tc_methanol-T)/(Tc_methanol-Tn_methanol))^0.38;
    Hv_ethanol = Hv_ethanol*((Tc_ethanol-T)/(Tc_ethanol-Tn_ethanol))^0.38;
    Hv = Hv_methanol*z+Hv_ethanol*(1-z);
    q = -CpV*(T-Tdp)/Hv;
end

sys = q;

% end mdlOutputs